#' ---
#' title: "Agenda Seeding: Spatial Panel Models"
#' date: "`r Sys.Date()`"
#' output: html_document
#' header-includes:
#'  - \usepackage{booktabs}
#'  - \usepackage{longtable}
#'  - \usepackage{array}
#'  - \usepackage{multirow}
#'  - \usepackage{wrapfig}
#'  - \usepackage{float}
#'  - \usepackage{colortbl}
#'  - \usepackage{pdflscape}
#'  - \usepackage{tabu}
#'  - \usepackage{threeparttable}
#'  - \usepackage{threeparttablex}
#'  - \usepackage[normalem]{ulem}
#'  - \usepackage{makecell}
#'  - \usepackage{dcolumn}
#'  - \usepackage{setspace}\doublespacing
#' ---



## ---- spatial_prep_spin_code, eval = FALSE, include = FALSE ----
# spin code to output Rmd / Rnw
# set output_format to "html_document" for html
# rmarkdown::render(input = here::here("code/spatial_panel_models3.R"), output_format = "pdf_document", clean = TRUE)




## ---- spatial_load_packages, include = FALSE  ----

library(tidyverse)
library(here)
library(naniar)
library(ggplot2)
library(stringr)
library(broom)

# spatial
library(maptools)
library(spdep)
library(rgdal)

# models
library(plm)
library(splm)

# tables
library(stargazer)


## ---- spatial_custom_functions, include = FALSE  ----

cache_logical <- FALSE

kable_format <- case_when(
    knitr::is_latex_output() ~ "latex",
    knitr::is_html_output()  ~ "html",
    TRUE                     ~ "markdown"
)

star_format <- case_when(
    knitr::is_latex_output() ~ "latex",
    knitr::is_html_output()  ~ "html",
    TRUE                     ~ "text"
)

# for stargazer
star.cut.vector <- c(0.05, NA, NA)


fiver <- function(x) {
    str_pad(
        string = x,
        width  = 5,
        side   = "left",
        pad    = "0"
    )
}



## ---- load_data_updated, eval = TRUE, include = FALSE  ----

load(here("data/voting_census_carter_rain_dca.Rdata"), verbose = TRUE)

# updated pdata
pdata_nv  <- vc2 %>% 
    rename(p_effect = dca_nv_bin) %>% 
    select(-matches('^car_|^dca_|^pwall|^prep|^pto|^pttl|^poth')) %>% 
    select(fips, year, everything()) %>% 
    na.omit()

pdata_v   <- vc2 %>% 
    rename(p_effect = dca_v_bin) %>% 
    select(-matches('^car_|^dca_|^pwall|^prep|^pto|^pttl|^poth')) %>% 
    select(fips, year, everything()) %>% 
    na.omit()

pdata_v2  <- vc2 %>% 
    rename(p_effect = car_bin) %>% 
    select(-matches('^car_|^dca_|^pwall|^prep|^pto|^pttl|^poth')) %>% 
    select(fips, year, everything()) %>% 
    na.omit()



## ---- load_formulas, include = FALSE  ----
form_bivariate <- pdem ~ p_effect

form_multiple  <- pdem ~ p_effect + locgov_exp_pc_log + per_hsplus + per_black + per_black_sq + med_age + med_income + per_unemp + per_urban + tot_pop_log + per_owner_housing + per_pop_growth + per_pop_foreign + plag 

form_lag       <- pdem ~ p_effect + locgov_exp_pc_log + per_hsplus + per_black + per_black_sq + med_age + med_income + per_unemp + per_urban + tot_pop_log + per_owner_housing + per_pop_growth + per_pop_foreign + plag 


## ---- run_plm_within, include = FALSE  ----

# run plm.calc chunk 
# index = fips, year
pout_nv_bi <- plm(formula = form_bivariate, data = pdata_nv, model = "within") 
pout_nv    <- plm(formula = form_multiple,  data = pdata_nv, model = "within") 
pout_v_bi  <- plm(formula = form_bivariate, data = pdata_v,  model = "within") 
pout_v     <- plm(formula = form_multiple,  data = pdata_v,  model = "within") 
pout_v2_bi <- plm(formula = form_bivariate, data = pdata_v2, model = "within") 
pout_v2    <- plm(formula = form_multiple,  data = pdata_v2, model = "within") 


## ---- plm_table_county_fe, results = 'asis', echo = FALSE  ----

# plm models table
stargazer(
    pout_nv_bi,
    pout_nv,
    pout_v_bi,
    pout_v,
    pout_v2_bi,
    pout_v2,
    column.labels    = c(
        "Nonviolent Protests (Olzak data)",
        "Violent Protests (Olzak data)",
        "Violent Protests (Carter data)"
    ) ,
    font.size        = "scriptsize",
    align            = TRUE,
    star.cutoffs     = star.cut.vector,
    notes.append     = FALSE,
    notes            = "*$p<0.05$",
    digits           = 2,
    dep.var.caption  = "DV: County-level Democratic Presidential Vote Share",
    dep.var.labels   = c(""),
    model.numbers    = TRUE,
    column.separate  = c(2, 2, 2) ,
    float.env        = "table",
    covariate.labels = c(
        "Protest `Treatment'",
        "log(PC Local Gov Exp)",
        "\\% HS+ Educ",
        "\\% Black",
        "(\\% Black)$^2$",
        "Median Age",
        "Median Income (000s)",
        "\\% Unemployment",
        "\\% Urban",
        "log(Population)",
        "\\% Owner Occ Housing",
        "\\% Pop Growth",
        "\\% Pop Foreign",
        "Lagged Dem Vote Share"
    ),
    omit.stat        = c("f", "ser"),
    title            = "Panel models, County Fixed Effects",
    type             = 'latex'
)



## ---- load_map_data, include = FALSE  ----

# http://clarkdatalabs.github.io/mapping_R/
# See County_2010Census_DP1 at:
# https://www2.census.gov/geo/tiger/TIGER2010DP1/

counties <- rgdal::readOGR(here("data/county_census/County_2010Census_DP1.shp"))



## ---- subset_map_data, cache = cache_logical, include = FALSE ----


matching_counties <- which(counties$GEOID10 %in% pdata_nv$fips)

length(matching_counties)


us48 <- counties[matching_counties, ]
dim(us48)

# poly to neighbor
nbq <- poly2nb(us48, queen=TRUE, row.names=us48$GEOID10)

# neighbor ids
nb_ids <- attr(nbq, "region.id")

length(nb_ids)

# convert nb object to listw object
usalw <- nb2listw(nbq)



## ---- create_panel, cache = cache_logical, include = FALSE ----
us48fips <- which(pdata_nv$fips %in% nb_ids)
p3 <- pdata_nv[us48fips, ] 

# dichotomize
p3$p_effect[p3$p_effect > 0] <- 1



## ---- spgm_bi, eval = TRUE, cache = cache_logical, include = FALSE  ----
sp_nv_bi <- spgm(
    formula       = form_bivariate,
    data          = p3,
    model         = "within",
    listw         = usalw,
    moments       = "initial",
    spatial.error = TRUE
)

#summary(sp_nv_bi)


## ---- spgm_multi, cache = cache_logical, include = FALSE ----
sp_nv <- spgm(
    formula       = form_multiple,
    data          = p3,
    model         = "within",
    listw         = usalw,
    moments       = "initial",
    spatial.error = TRUE
)

#summary(sp_nv)



## ---- create_panel_v, cache = cache_logical, include = FALSE ---- 
us48fips <- which(pdata_v$fips %in% nb_ids)
p3       <- pdata_v[us48fips, ]

# dichotomize
p3$p_effect[p3$p_effect > 0] <- 1



## ---- spgm_bi_v, eval = TRUE, cache = cache_logical, include = FALSE ---- 
sp_v_bi <- spgm(
    formula       = form_bivariate,
    data          = p3,
    model         = "within",
    listw         = usalw,
    moments       = "initial",
    spatial.error = TRUE
) 

#summary(sp_v_bi)


## ---- spgm_multi_v, cache = cache_logical, include = FALSE ---- 
sp_v <- spgm(
    formula       = form_multiple,
    data          = p3,
    model         = "within",
    listw         = usalw,
    moments       = "initial",
    spatial.error = TRUE
) 

#summary(sp_v)



## ---- create_panel_v2, cache = cache_logical, include = FALSE ---- 
us48fips <- which(pdata_v2$fips %in% nb_ids)
p3       <- pdata_v2[us48fips, ]

# dichotomize
p3$p_effect[p3$p_effect > 0] <- 1



## ---- spgm_bi_v2, eval = TRUE, cache = cache_logical, include = FALSE ---- 
sp_v2_bi <- spgm(
    formula       = form_bivariate,
    data          = p3,
    model         = "within",
    listw         = usalw,
    moments       = "initial",
    spatial.error = TRUE
)

#summary(sp_v2_bi)


## ---- spgm_multi_v2, cache = cache_logical, include = FALSE ---- 
sp_v2 <- spgm(
    formula       = form_multiple,
    data          = p3,
    model         = "within",
    listw         = usalw,
    moments       = "initial",
    spatial.error = TRUE
) 

#summary(sp_nv)



## ---- extract_model_stats_to_pass_to_stargazer, include = FALSE ----

# stargazer can't take spgm models so building by hand

N     <- c(length(sp_nv$residuals), length(sp_v$residuals), length(sp_v2$residuals))
rho   <- c(sp_nv$rho[1], sp_v$rho[1], sp_v2$rho[1] )
sigma <- c(sp_nv$rho[2], sp_v$rho[2], sp_v2$rho[2] )




## ---- create_table, results = 'asis' ---- 
# Bivariate splm models
# prints Protest treatment in Constant row??
stargazer(pout_nv_bi, pout_v_bi, pout_v2_bi,
          coef = list(
              summary(sp_nv_bi)$CoefTable[, 1],
              summary(sp_v_bi )$CoefTable[, 1],
              summary(sp_v2_bi)$CoefTable[, 1]
          ),
          se = list(
              summary(sp_nv_bi)$CoefTable[, 2],
              summary(sp_v_bi )$CoefTable[, 2],
              summary(sp_v2_bi)$CoefTable[, 2]
          ),
          font.size        = "scriptsize",
          omit.stat        = c("all"),
          #omit             = c("Constant"),
          add.lines        = list(c("N", N) ,
                                 c("$\\rho$", round(rho, 2)) ,
                                 c("$\\sigma^2_v$", round(sigma, 3))),
          align            = TRUE,
          star.cutoffs     = star.cut.vector,
          notes.append     = FALSE,
          notes            = "*$p<0.05$",
          digits           = 2,
          dep.var.caption  = "DV: County-level Change in Democratic Vote Share",
          dep.var.labels   = c(""),
          model.numbers    = TRUE,
          column.labels    = c(
              "Nonviolent Protests (Olzak data)",
              "Violent Protests (Olzak data)",
              "Violent Protests (Carter data)"
          ),
          column.separate  = c(1, 1, 1) ,
          float.env        = "table",
          covariate.labels = c("Protest `Treatment'"),
          title            = "Spatial Panel Models (Bivariate)",
          type             = 'latex')



## ---- final_table_save, results = 'asis' ---- 
# Multivariate splm models
#spgm_multi_table_out <- 
# write to file with capture.output
capture.output( 
    stargazer(
    pout_nv,
    pout_v,
    pout_v2,
    coef = list(
        summary(sp_nv)$CoefTable[, 1],
        summary(sp_v )$CoefTable[, 1],
        summary(sp_v2)$CoefTable[, 1]
    ),
    se = list(
        summary(sp_nv)$CoefTable[, 2],
        summary(sp_v )$CoefTable[, 2],
        summary(sp_v2)$CoefTable[, 2]
    ),
    font.size        = "scriptsize",
    omit.stat        = c("all"),
    omit             = c("Constant"),
    add.lines        = list(c("N", N) ,
                            c("$\\rho$", round(rho, 2)) ,
                            c("$\\sigma^2_v$", round(sigma, 2))),
    align            = TRUE,
    star.cutoffs     = star.cut.vector,
    notes.append     = FALSE,
    notes            = "*$p<0.05$",
    digits           = 2,
    dep.var.caption  = "DV: County-level Change in Democratic Vote Share",
    dep.var.labels   = c(""),
    model.numbers    = TRUE,
    column.labels    = c(
        "Nonviolent Protests (Olzak data)",
        "Violent Protests (Olzak data)",
        "Violent Protests (Carter data)"
    ),
    column.separate  = c(1, 1, 1) ,
    float.env        = "table",
    covariate.labels = c(
        "Protest `Treatment'",
        "log(PC Local Gov Exp)",
        "\\% HS+ Educ",
        "\\% Black",
        "(\\% Black)$^2$",
        "Median Age",
        "Median Income (000s)",
        "\\% Unemployment",
        "\\% Urban",
        "log(Population)",
        "\\% Owner Occ Housing",
        "\\% Pop Growth",
        "\\% Pop Foreign",
        "Lagged Dem Vote Share"
    ),
    title = "Spatial Panel Models",
    type  = "latex"
),
    file = here("code/spatial_panel_multi_table.tex") )




## ---- save_spatial_regression_results, include = FALSE ----
#

sp_coef <- c(summary(sp_nv)$CoefTable[1, 1],
             summary(sp_v )$CoefTable[1, 1],
             summary(sp_v2)$CoefTable[1, 1])
sp_se   <- c(summary(sp_nv)$CoefTable[1, 2],
             summary(sp_v )$CoefTable[1, 2],
             summary(sp_v2)$CoefTable[1, 2])

sp_coef
sp_se

save(pout_nv, pout_v, pout_v2, sp_nv, sp_v, sp_v2, sp_coef, sp_se, 
     file = here("data/spatial_panel_multi2.Rdata"))


## ---- final_table_in_doc, results = 'asis' ---- 

stargazer(
    pout_nv,
    pout_v,
    pout_v2,
    coef = list(
        summary(sp_nv)$CoefTable[, 1],
        summary(sp_v )$CoefTable[, 1],
        summary(sp_v2)$CoefTable[, 1]
    ),
    se = list(
        summary(sp_nv)$CoefTable[, 2],
        summary(sp_v )$CoefTable[, 2],
        summary(sp_v2)$CoefTable[, 2]
    ),
    font.size        = "scriptsize",
    omit.stat        = c("all"),
    omit             = c("Constant"),
    add.lines        = list(c("N", N) ,
                            c("$\\rho$", round(rho, 2)) ,
                            c("$\\sigma^2_v$", round(sigma, 2))),
    align            = TRUE,
    star.cutoffs     = star.cut.vector,
    notes.append     = FALSE,
    notes            = "*$p<0.05$",
    digits           = 2,
    dep.var.caption  = "DV: County-level Change in Democratic Vote Share",
    dep.var.labels   = c(""),
    model.numbers    = TRUE,
    column.labels    = c(
        "Nonviolent Protests (Olzak data)",
        "Violent Protests (Olzak data)",
        "Violent Protests (Carter data)"
    ),
    column.separate  = c(1, 1, 1) ,
    float.env        = "table",
    covariate.labels = c(
        "Protest `Treatment'",
        "log(PC Local Gov Exp)",
        "\\% HS+ Educ",
        "\\% Black",
        "(\\% Black)$^2$",
        "Median Age",
        "Median Income (000s)",
        "\\% Unemployment",
        "\\% Urban",
        "log(Population)",
        "\\% Owner Occ Housing",
        "\\% Pop Growth",
        "\\% Pop Foreign",
        "Lagged Dem Vote Share"
    ),
    title = "Spatial Panel Models",
    #out   = "spatial_panel_multivariate_table.tex",
    type  = "latex"
)

